6 Georreferencia
Análisis Exploratorio de Datos | Licenciatura en Estadística | FCEyE | UNR
6.1 Introducción
La Georreferencia es una disciplina ligada al análisis de datos recolectados a nivel geográfico. Esta característica implica que la información que poseemos puede visualizarse de manera más eficiente desde una perspectiva espacial, es decir, representando los datos en mapas y no mediante gráficos, tablas u otros instrumentos.
Ejemplo: deseamos representar los casos confirmados de Covid-19 en Argentina. Existen dos enfoques desde los cuales se pueden analizar estos datos geográficamente:
Mapas Coropléticos (choropleth maps): utilizan paletas de colores para representar una variable a nivel geográfico, asignando un color a cada unidad representada (país, provincia, departamento, etc.) en base a los valores que le corresponden. Las variables de interés en estos casos pueden ser de tipo cualitativo (por ej. nivel de riesgo epidemiológico) o bien cuantitativo (cantidad de casos, tasa de positividad en tests de Covid-19, etc.). Bajo este enfoque, lo geográfico se impone únicamente como medio de visualización y no como variable de interés primario.
Marcadores: consiste en señalar coordenadas específicas mediante un punto o marcador, para mostrar dónde ocurren los eventos que estamos estudiando. En estos casos el interés está puesto en visualizar qué tan comunes son ciertos eventos y cómo se ubican espacialmente. En el ejemplo de Covid-19, podemos marcar los domicilios de las personas contagiadas para encontrar focos de infección a partir de la distribución geográfica de los puntos.
En esta clase repasaremos casos de aplicación de ambos enfoques.
6.1.1 Proyecciones
Un concepto fundamental a la hora de trabajar con este tipo de datos es el de proyección geográfica. Todo punto que se encuentre sobre un plano puede representarse mediante coordenadas, consistentes en un par de números \((x,y)\) que miden la distancia de ese punto con respecto al origen.
La forma en que se calculan las distancias, y la ubicación del origen, dependen del Sistema de Referencia de Coordenadas (CRS) utilizado. El CRS más común es WGS84, el cual representa a la Tierra como un elipsoide y utiliza las siguientes definiciones:
- Latitud: coordenada que especifica la ubicación norte-sur de cualquier punto sobre la tierra. Varía desde 90° en el polo norte hasta -90° en el polo sur. A las líneas de latitud constante se las llama “paralelos”. La latitud 0° corresponde a línea del Ecuador, la cual divide al planeta en los hemisferios norte y sur.
- Longitud: coordenada que especifica la ubicación este-oeste de cualquier punto sobre la tierra. A las líneas de longitud constante se las llama “meridianos”. La longitud 0° corresponde al meridiano de Greenwich, en Inglaterra; los restantes meridianos se miden hasta el 180° (dirección este) y -180° (dirección oeste), uniéndose sobre la Línea Internacional del Cambio de Fecha (anti-meridiano), ubicada sobre el océano Pacífico.
Dado que la tierra es esférica, representar latitudes y longitudes usando ejes cartesianos no es recomendable, ya que distorsiona mucho las áreas, formas y distancias. Para sobrellevar este problema se aplican proyecciones no lineales que tratan de minimizar estas distorsiones (no es posible conservar todas las propiedades al mismo tiempo; por ejemplo, si respetamos las distancias reales, las formas y las superficies de las regiones se distorsionan).
Debajo se muestran algunos ejemplos de proyecciones:
- Una de las proyecciones más populares es la de Mercator. Si bien fue utilizada por Google Maps hasta 2018, posee grandes distorsiones en cuanto al área de las regiones representadas. Por ejemplo, vemos a Groenlandia aproximadamente del mismo tamaño que América del Sur, cuando en realidad sus superficies son de aproximadamente 2 y 18 millones de \(km^2\) respectivamente:
6.1.2 Georreferencia con R
En los últimos años se produjo un boom en el uso de Sistemas de Información Geográfica (SIGs) para el análisis de información espacial. Si bien programas como ArcGIS o QGIS son muy populares en este contexto, autores como Brunsdon & Comber (2019) y Lovelace et al. (2019) remarcan que la reproducibilidad de los análisis se ve comprometida debido a las interfaces gráficas utilizadas por esos productos.
Puntualmente, las interfaces amigables para el usuario (user-friendly), con sus menús y opciones clickeables, atentan contra el trabajo colaborativo y la replicación del proceso de investigación. Por este motivo creemos que es importante trabajar estos datos desde la perspectiva del lenguaje R.
El ecosistema de R ha sido un terreno extremadamente fértil para el desarrollo y crecimiento de paquetes ligados al análisis geocomputacional. En este link se listan unos 200 paquetes referidos a esta temática, entre los que se destacan
sp,gstat,rgdal,stars,maptools,spatial,raster,mapviewyRgoogleMaps, entre tantos otros.En particular, nosotros emplearemos los siguientes paquetes a lo largo de la clase:
#Mapas
library(sf)
library(tmap)
library(leaflet)
library(spData)
#Lectura y manipulación de datos
library(readxl)
library(readr)
library(dplyr)
library(tidyr)
library(stringr)
library(ggplot2)
library(htmltools)Este apunte está organizado de la siguiente manera:
- Comenzaremos nuestro recorrido repasado los recursos disponibles online para descarga de mapas.
- Luego nos centraremos en el paquete
sf, el cual nos brinda la posibilidad de manipular datos geoespaciales, integrándolo conggplot2para producir mapas estáticos. - Por último, estudiaremos las funciones más importantes de
tmapyleaflet, capaces de generar mapas dinámicos que le otorgan un salto de calidad a nuestro trabajo.
6.2 Buscando Mapas en la Web
- Existen muchas páginas desde donde podemos descargar archivos con información geoespacial para utilizar en nuestros mapas. A continuación repasaremos algunas de ellas.
6.2.1 INDEC
En la web del Instituto Nacional de Estadística y Censos (INDEC) de la República Argentina hay una sección dedicada a la Cartografía. Desde allí podemos descargar información con distintos niveles de jerarquía o subdivisiones:
- Mapas provinciales
- Mapas departamentales
- Localidades censales
- Radios censales
- Etc.
- En general, las divisiones geográficas empleadas por INDEC resultan de interés a la hora de analizar datos recolectados durante los Censos de población, o bien datos pertenecientes a la Encuesta Permanente de Hogares (EPH).
6.2.2 IGN
El Instituto Geográfico Nacional que depende del Ministerio de Defensa de la República Argentina es la fuente de datos oficial en todo lo que respecta a la cartografía de nuestro país.
Su página web ofrece una gran variedad de recursos, entre los que se incluyen mapas de todo tipo, imágenes satelitales, capas de información geoespacial, etc.
El archivo de formato shapefile que usaremos en los ejemplos presentados a continuación puede descargarse desde la web del IGN. Posee un peso de aproximadamente 34 MB y contiene los límites de todas las subdivisiones administrativas (departamentos o partidos) de las provincias argentinas.
6.2.3 GADM
Una web que ofrece mapas políticos de casi todos los países y territorios a nivel mundial es GADM. Si bien es un recurso gratuito, sus mapas poseen una resolución mayor a la de software pagos como ArcGIS.
Existen tres niveles de mapas en GADM:
- Países
- Países con provincias
- Países con sus provincias y subdivisiones menores (departamentos, partidos, etc.)
6.3 Mapas Estáticos con sf
Como siempre, el primer paso en cualquier proceso de análisis de datos consiste en importar la base. Para el caso de información geográfica, podemos emplear la función
read_sf()del paquetesf.Recordemos que el archivo a leer posee formato .shp (shapefile) y contiene los límites de cada departamento (subdivisiones provinciales) de nuestro país. Suponiendo que este archivo se llama departamento y se encuentra almacenado en una carpeta llamada argentina dentro de nuestro directorio de trabajo, la sintaxis a emplear para importarlo es la siguiente:
mi_mapa <- read_sf(
dsn = "../data/unidad03/argentina", #dsn: carpeta donde está el archivo shp
layer = "departamento" #layer: nombre del archivo shp
)- Para los ejemplos que vienen a continuación, trabajaremos únicamente con los departamentos de la Provincia de Santa Fe:
#fdc: nombre del ente gubernamental que recolectó la información geográfica
#cada provincia posee uno propio
santafe <- mi_mapa %>%
filter(fdc == "Servicio de Catastro e Información Territorial") %>%
mutate(Nombre = str_wrap(nam, 8)) #separamos nombres largos en 2 lineas- Analicemos la estructura de este objeto:
str(santafe)sf [19 × 10] (S3: sf/tbl_df/tbl/data.frame)
$ gid : int [1:19] 175 176 170 171 172 173 174 429 430 431 ...
$ objeto : chr [1:19] "Departamento" "Departamento" "Departamento" "Departamento" ...
$ fna : chr [1:19] "Departamento Vera" "Departamento 9 de Julio" "Departamento Las Colonias" "Departamento San Justo" ...
$ gna : chr [1:19] "Departamento" "Departamento" "Departamento" "Departamento" ...
$ nam : chr [1:19] "Vera" "9 de Julio" "Las Colonias" "San Justo" ...
$ in1 : chr [1:19] "82133" "82077" "82070" "82112" ...
$ fdc : chr [1:19] "Servicio de Catastro e Información Territorial" "Servicio de Catastro e Información Territorial" "Servicio de Catastro e Información Territorial" "Servicio de Catastro e Información Territorial" ...
$ sag : chr [1:19] "IGN" "IGN" "IGN" "IGN" ...
$ geometry:sfc_MULTIPOLYGON of length 19; first list element: List of 1
..$ :List of 1
.. ..$ : num [1:75079, 1:2] -60.2 -59.8 -59.8 -59.8 -59.8 ...
..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
$ Nombre : chr [1:19] "Vera" "9 de\nJulio" "Las\nColonias" "San\nJusto" ...
- attr(*, "sf_column")= chr "geometry"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA
..- attr(*, "names")= chr [1:9] "gid" "objeto" "fna" "gna" ...La última columna del conjunto de datos, llamada geometry, es un objeto de clase
sfc_MULTIPOLYGONdonde se guardan las coordenadas que representan los límites de cada región.La primera fila del dataset, correspondiente al departamento de Vera, está conformada por una lista de 75.079 pares de puntos \((x,y)\) que, al ser conectados, construyen el contorno del departamento.
Ejemplo para el departamento Rosario:
rosario <- santafe$geometry[[11]] %>% as.matrix() %>% as.data.frame()
head(rosario) V1 V2
1 -60.80908 -32.82947
2 -60.80886 -32.83369
3 -60.80596 -32.83359
4 -60.79912 -32.83334
5 -60.79633 -32.83323
6 -60.79540 -32.83319ggplot(data = rosario) +
aes(x = V1, y = V2) +
geom_point()
Para poder graficar estos datos de manera directa, sin necesidad de extraer las coordenadas, podemos usar en conjunto los paquetes
ggplot2ysf.Una de las tantas capas que podemos agregar a un
ggplotesgeom_sf(), la cual toma un objeto con información espacial (en este caso la columna con clasesfc_MULTIPOLYGON) y une los puntos mediante líneas, para graficar los límites geográficos de cada región.Veamos a cada uno de los departamentos de Santa Fe representados por un color diferente, asociado a su propio nombre (
fill = Nombre):
#No hace falta mencionar a la columna 'geometry' ya que es el
#nombre estándar donde se guarda la información geográfica en este tipo de bases
ggplot(data = santafe) +
aes(fill = nam) +
geom_sf() +
geom_sf_text(aes(label = Nombre), size = 2) +
theme_bw() +
theme(legend.position = "none")
6.3.1 Datos de cultivos
En la visualización generada arriba la variable representada es categórica (nombre del departamento). Otra posible manera de crear mapas coropléticos es asignar colores a variables numéricas, por ejemplo algunas de las que se encuentran en la base de cultivos que analizamos en clases anteriores.
¿Cómo podemos implementar esta idea? En primer lugar, importamos y manipulamos el conjunto de datos para elegir sólo algunas variables de interés:
cultivos <- read_excel("../data/unidad03/cultivos.xlsx", na = "-")
soja_sf <- cultivos %>%
filter(
campaña == "2019/20",
prov == "SANTA FE",
cultivo == "Soja total"
) %>%
mutate(Soja = prod/1000) %>%
select(Departamento = dpto, Soja)- La base resultante posee 19 filas, una para cada departamento santafesino, y datos sobre la producción de soja (en miles de toneladas) durante la campaña 2019/20:
| Departamento | Soja |
|---|---|
| 9 DE JULIO | 269.220 |
| BELGRANO | 530.474 |
| CASEROS | 759.772 |
| CASTELLANOS | 1030.429 |
| CONSTITUCION | 561.055 |
| GARAY | 9.000 |
| GENERAL LOPEZ | 2024.850 |
| GENERAL OBLIGADO | 86.490 |
| IRIONDO | 677.818 |
| LA CAPITAL | 81.200 |
| LAS COLONIAS | 667.804 |
| ROSARIO | 196.280 |
| SAN CRISTOBAL | 414.891 |
| SAN JAVIER | 20.160 |
| SAN JERONIMO | 449.692 |
| SAN JUSTO | 265.702 |
| SAN LORENZO | 323.504 |
| SAN MARTIN | 982.996 |
| VERA | 48.630 |
- El segundo paso consiste en unir esta información con aquella presente en el shapefile, lo cual puede hacerse usando la función
left_join()dedplyr. Dado que la llave (key) para unificar bases es el nombre del departamento, previamente llevamos a cabo un proceso de unificación de escrituras:
#Comparacion de nombres entre ambas bases:
#cbind(sort(santafe$nam), sort(soja_sf$Departamento))
sf_unido <- santafe %>%
mutate(
Departamento = str_to_upper(nam), #pasamos todos los nombres a mayúscula
Departamento = chartr("ÁÉÍÓÚ", "AEIOU", Departamento) #borramos tildes
) %>%
left_join(soja_sf)- Ya estamos en condiciones de comparar producciones de cultivos para los departamentos santafesinos. Veamos el caso de la Soja:
ggplot(data = sf_unido) +
aes(fill = Soja) +
geom_sf() +
geom_sf_text(aes(label = str_wrap(nam, 8)), size = 2) +
scale_fill_gradient(low = "#E0ECF8", high = "#0174DF") +
theme_bw()
- Cualquier modificación que deseemos aplicar sobre un mapa como este puede llevarse a cabo siguiendo las reglas correspondientes a gráficos comunes y corrientes producidos por
ggplot2.
6.3.2 Trabajo en Equipo
- Para practicar la creación de mapas vamos a usar el conjunto de datos world del paquete
spData. Cada una de las 177 filas de esta base corresponde a un país diferente, para los cuales se registran diversas variables: nombre, continente, superficie, población, etc.
data(world) #Cargamos la base
head(world) #Vemos sus primeras 6 filasSimple feature collection with 6 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
Geodetic CRS: WGS 84
# A tibble: 6 × 11
iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp
<chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
1 FJ Fiji Oceania Oceania Melanesia Sove… 1.93e4 8.86e5 70.0
2 TZ Tanzania Africa Africa Eastern … Sove… 9.33e5 5.22e7 64.2
3 EH Western S… Africa Africa Northern… Inde… 9.63e4 NA NA
4 CA Canada North Am… Americas Northern… Sove… 1.00e7 3.55e7 82.0
5 US United St… North Am… Americas Northern… Coun… 9.51e6 3.19e8 78.8
6 KZ Kazakhstan Asia Asia Central … Sove… 2.73e6 1.73e7 71.6
# ℹ 2 more variables: gdpPercap <dbl>, geom <MULTIPOLYGON [°]>
Usando funciones del paquete
sf, graficar los países de América del Sur asignando colores de acuerdo a su esperanza de vida. La paleta debe asignar el color rojo a valores bajos y el verde a valores altos de esta variable.Resultado esperado:
6.4 Mapas Dinámicos
Los mapas estáticos, como los que estuvimos viendo hasta ahora, son útiles pero en cierto sentido limitados: sólo nos muestran una imagen fija, para un período de tiempo dado, y no nos permiten cambiar la información graficada.
Por suerte para nosotros, R ofrece un variado menú de opciones a la hora de dar vida a nuestros mapas, permitiéndonos hacerlos dinámicos o interactivos. Esto implica que podemos modificar las variables o el área a visualizar, tal como lo hacemos habitualmente con Google Maps o servicios similares.
6.4.1 Paquete tmap
Una de las opciones para crear este tipo de mapas en R es a través del paquete
tmap. Esta librería emplea una sintaxis similar a la deggplot2. Existe una analogía entre las funciones de ambos paquetes; por ejemplo,ggplot()se convierte entm_shape()para dar inicio al gráfico, las capas de tipogeom_se llaman ahoratm_, etc.Veamos un ejemplo simple donde aplicamos
tm_shape()para obtener un mapa coroplético equivalente al anterior, donde comparamos la producción de soja para diferentes departamentos de Santa Fe.Algunos detalles sobre el código:
- Usamos
tmap_mode("view")para activar la funcionalidad interactiva de esta librería; si deseamos mapas estáticos ejecutamostmap_mode("plot"). - La variable agregada con nombre ID almacena el texto que se mostrará en el mapa cada vez que pasamos el mouse sobre una región determinada (se le puede dar formato usando lenguaje HTML).
- La función
tm_basemap()se usa para elegir el tipo de mapa de fondo utilizado. La lista completa de opciones puede consultarse en esta web.
- Usamos
tmap_mode("view")
sf_unido %>%
mutate(ID = paste0("Depto. ", nam, " - ", Soja, " millones de tn.")) %>%
tm_shape() +
tm_fill(id = "ID", col = "Soja", palette = "Oranges", alpha = 0.75) +
tm_borders("red", lwd = 2) +
tm_minimap() +
tm_basemap("OpenStreetMap")